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Abstract 

Monte Carlo is a simple and flexible tool that is widely used in 
computational finance. In this context, it is common for the quantity 
of interest to be the expected value of a random variable defined via 
a stochastic differential equation. In 2008, Giles proposed a remark¬ 
able improvement to the approach of discretizing with a numerical 
method and applying standard Monte Carlo. His multilevel Monte 
Carlo method offers a speed up of 0(e -1 ), where e is the required 
accuracy. So computations can run 100 times more quickly when two 
digits of accuracy are required. The “multilevel philosophy” has since 
been adopted by a range of researchers and a wealth of practically sig¬ 
nificant results has arisen, most of which have yet to make their way 
into the expository literature. In this work, we give a brief, accessible, 
introduction to multilevel Monte Carlo and summarize recent results 
applicable to the task of option evaluation. 

Keywords computational complexity, control variate, Euler-Maruyama, 
Monte Carlo, option value, stochastic differential equation, variance reduc¬ 
tion. 


1 Aims 

Finding the appropriate market value of a financial option can usually be 
formulated as an expected value computation pO] !23 ; , [38]. In the case where 
the product underlying the option is modelled as a stochastic differential 
equation (SDE), we may 
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• simulate the SDE numerically to compute many independent sample 
paths, and then 

• combine the option payoff from each path in order to obtain a Monte 
Carlo estimate, and an accompanying confidence interval. 

Compared with other approaches, notably the direct discretization of a par¬ 
tial differential equation based formulation of the problem, a Monte Carlo 
computation has the advantages of (a) being simple to implement and (b) 
being flexible enough to cope with a wide range of underlying SDE models 
and option payoffs. On the downside, Monte Carlo is typically expensive in 
terms of computation time 2(). [23]. 

In the seminal 2008 paper hi, Giles pulled together ideas from numer¬ 
ical analysis, stochastic analysis and applied statistics in order to deliver a 
dramatic improvement on the efficiency of the “SDE simulation plus Monte 
Carlo” approach. If the required level of accuracy, in terms of confidence 
interval, is e, the multilevel approach essentially improves the computational 
complexity by a factor of e. So for a calculation requiring two digits of 
accuracy, we obtain a hundredfold improvement in computation time. Mul¬ 
tilevel Monte Carlo has rapidly become an extremely hot topic in the field of 
stochastic computation, impacting on a wide range of application areas. In 
particular, technical reviews of research progress in the field have begun to 
appear [16, T8j and a comprehensive survey is currently in progress by Giles 
for the journal Acta Numerica. However, the area is still sufficiently new 
that most textbooks in computational finance do not introduce the topic, 
and hence it has not been fully integrated into typical graduate-level classes 
and development courses for practitioners. For this reason, we present here 
an accessible introduction to the multilevel Monte Carlo approach, and give 
a brief overview of the current state of the art with respect to financial option 
valuation. 

In section [2] we discuss the underlying SDE simulation. Section [3] then 
considers the complexity of standard Monte Carlo in this setting. In section [4] 
we give some motivation for the multilevel approach, which is introduced and 
analysed in section [5j Section [6] illustrates the performance of the algorithm 
in practice, using code that has been made available by Giles. In section [7] 
we give pointers to multilevel research in option valuation that has built on 
[2]. Section [8] concludes with a brief discussion. 
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2 Convergence in SDE Simulation 

We consider an Ito SDE of the form 


dX{t) = f(X(t))dt + g(X(t))dW (£), X(0) = X 0 . (1) 

Here, / : M m —> W 1 and g : M m — > M mxd are given functions, known as the 
drift and diffusion coefficients, respectively, and W ( t ) e is standard Brow¬ 
nian motion. The initial condition Xq is supplied and we wish to simulate 
the SDE over the fixed time interval [0,T]. The Euler-Maruyama method 
pH! [53] computes approximations X n « X(t n ), where t n = nAt, according 
to X 0 = A(0) and, for n — 1,2,... N — 1, 

X n+1 = X n + f(X n )At + g(X n )AW n , (2) 

where At = T/N is the stepsize and A W n = W{t n+ 1 ) — W(t n ) is the relevant 
Brownian motion increment. 

In the study of the accuracy of SDE simulation methods, the two most 
widely used convergence concepts are referred to as weak and strong [31], [33]. 
Roughly, 

• weak convergence controls the error of the means, whereas, 

• strong convergence controls the mean of the error. 

To prove weak and strong convergence results, we must impose conditions 
on the SDE. For example it is standard to assume that / and g in (JTj) satisfy 
global Lipschitz conditions; that is, there exists a constant L such that 


l/W-/(y)l < L\x~v\ and \g(x) - g(y)\ < L\x - y\, for all x, y 6 K m . 

(3) 

Here and throughout we take || ■ || to be the Euclidean norm. Under such con¬ 
ditions, and for appropriate initial data, it follows that the Euler-Maruyama 
method has weak order one, so that 


sup (E[X(t n )]-K[X n ]) = 0(At). (4) 

0 <t n <T 

In the sense of strong error, which involves the mean of the absolute differ¬ 
ence between the two random variables at each grid point, Euler-Maruyama 
achieves only an order of one half in general: 


E 


sup | X(t n )-X n 
g<t„<T 


O(AU). 


(5) 
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More generally, for any m > 1 and sufficiently small At there is a constant 
C = C(m) such that 


E 


sup | X(t n )-X n 

0 <t n <T 


< CAt m/2 . 


( 6 ) 


Strong convergence is sometimes described as a pathwise property. This can 
be understood via the Borel-Cantclli Lemma. For example, in [3Uj it is shown 
that given any e > 0 there exists a path-dependent constant K = K (e) such 
that, for all sufficiently small At, 

sup \X(t n ) - X n \ < K(e)h^~ e . 

0 <t n <T 

In the setting of this work, it is tempting to argue that strong convergence 
is not relevant; if we wish to compute an expected value based on the SDE 
solution then following individual paths accurately is not important. How¬ 
ever, we will see in section [5] that the analysis in [ITj justifying multilevel 
Monte Carlo makes use of both weak and strong convergence properties. 

To conclude this section, we remark that the analysis of SDE simulation 
on problems that violate the global Lipschitz conditions ([3]) is far from com¬ 
plete. In the case of SDE models for financial assets and interest rates, issues 
may arise through faster than linear growth at infinity and also through un¬ 
bounded derivatives at the origin. For example, both complications occur in 
the class of scalar interest rate models from pQ, 

dX{t) = [a^Xity 1 - a 0 + a x X{t) - a 2 X(t ) r ) dt + a 3 X(t) p dW(t), 

where the a, are positive constants and r, p > 1. Although some positive 
results are available for specific nonlinear structures [23, ESI 126’ [39], there 
has also been a sequence of negative results showing how Euler-Maruyama 
can break down on nonlinear SDEs [23, 2T; [33]. 


3 SDE Simulation and Standard Monte Carlo 

Given the SDE (JT|) , suppose we wish to approximate the final time expected 
value of the solution, E[X(T)], using Monte Carlo with Euler-Maruyama. 
We will let e denote the accuracy requirement in terms of confidence interval 
width; fixing on a 95% confidence level to be concrete, we therefore wish to 
be in a position where applying the algorithm independently a large number 
of times, the exact answer would be within ±e of our computed answer with 
frequency at least 0.95. 
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Let X \y denote the Euler-Maruyama final time approximation along the 
sth path. Using M Monte Carlo samples we may form the sample average 

1 


M 


a M ~ TT X 


M 


S=1 


The overall error in onr approximation has the form 

a M -E[X(T)) = a M -E[X(T)-X N + X N ] 

= a M ~ E[A*] + E[X N - X(T)}. 


(7) 


Note that Xn denotes a random variable describing the result of applying 
Euler-Maruyama §, whereas each xjy is an independent sample from the 
distribution given by Xn- The expression ([7]) breaks down the error into 
two terms. The statistical error, aM — E [X n ], is concerned with how well 
we can approximate an expected value from a finite number of samples; 
it does not depend on how accurately the numerical method approximates 
the SDE (in particular it does not depend significantly on At) and it will 
generally decrease if we take more sample paths. The discretization error, or 
bias, E[X/v — X(T)], arises because we have approximated the SDE with a 
difference equation; this is the discrepancy that would remain if we had access 
to the exact expected value of the numerical solution and it will generally 
decrease if we reduce the stepsize. 

Standard results 1201 [37] tell us that the statistical error aM — E[X n ] 
can be described via a confidence interval of width 0(l/\/M). The weak 
convergence property Q shows that the bias E[X V — X(T )] behaves like 
O(Af); so we must add this amount to the confidence interval width. We 
arrive at an overall confidence interval of width 0(1/a/M) + O(At). To 
achieve our required target accuracy of e, we see that 1 /\[M and At should 
scale like e. In other words, M should scale like e -2 and At should scale like 
e. 

It is reasonable to measure computational cost by counting either the 
number of times that the drift and diffusion coefficients, / and g , are evalu¬ 
ated, or the number of times that a random number generator is called. In 
either case, the cost per path is proportional to 1/At, and hence the total 
cost of the computation scales like Mf At. We argued above that M should 
scale like e~ 2 and At should scale like e. Here is the conclusion: 

we may achieve accuracy e by combining Euler-Maruyama and 
standard Monte Carlo at an overall cost that scales like e~ 3 . 

One approach to improving the computational complexity is to replace 
Euler-Maruyama with a simulation method of higher weak order [U EU [33]. 
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If we use a second order method, so that Q is replaced by 
sup (E [X(t n )} - E[X n ]) = 0(Af 2 ), 

0 <t n <T 

then a straightforward adaption of the arguments above lead to the following 
conclusion: 

we may achieve accuracy e by combining a second order weak 
method and standard Monte Carlo at an overall cost that scales 
like e" 2 ' 5 . 

We note, however, that establishing second order weak convergence requires 
extra smoothness assumptions to be placed on the SDE coefficients. 

As we show in section[5j the method of Giles p3] has the following feature: 

we may achieve accuracy e by using Euler-Maruyama in a mul¬ 
tilevel Monte Carlo setting at an overall cost that scales like 
e- 2 (loge) 2 . 

Moreover, by using a higher strong order method, such as Milstein PU [33], 
it is possible to reduce the multilevel Monte Carlo cost to the order of e -2 

[Sl¬ 
it is worth pausing to admire an 0(e~ 2 ) computational complexity count. 
Suppose we are given an exact expression for the SDE solution, as a function 
of W{t). Hence, we are able to compute exact samples, without the need 
to apply a numerical method. A standard Monte Carlo approach requires 
1/y/M to scale like e in order to achieve the required confidence interval 
width. If we regard the evaluation of each exact X (T) sample as having 
0(1) cost, then the cost overall will be proportional to M; that is, e~ 2 . 
In this sense, with a multilevel approach the numerical analysis comes for 
free ; we can solve the problem as quickly as one for which we have an exact 
pathwise expression for the SDE solution. 


4 Motivation for the Multilevel Approach 


We can motivate the multilevel approach by considering a series expansion 
of Brownian motion, where the coefficients are random variables. The Paley- 
Wiener representation over [0, 2n\ has the form 


W(t) 



+ 


2 z sin (l nt ) 
n 


( 8 ) 
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where the {Zj}j> 0 are i.i.d. and N( 0,1); see, for example, [32]• In Figure[l]we 
draw samples for the Z % and plot the curves arising when the inhnite series 
in ([I]) is truncated to Yln=v f° r ^ = I) 2, 5, 10, 50 and 200. It is clear that 
the early terms in the series affect the overall shape, while the later terms 
add fine detail. From this perspective, it is intuitively reasonable that we 
can build up information at different resolution scales, with the finer scales 
having less impact on the overall picture. 

Now, we may view Monte Carlo as requiring a “black box” that returns 
independent samples. In our numerical SDE context, the samples come from 
a distribution that is only approximately correct, and the black box (the 
Euler-Maruyama method) comes with a dial. Turning the dial corresponds 
to changing At. Samples with a smaller At are more expensive—we have to 
wait longer for them because the paths contain more steps. The multilevel 
Monte Carlo algorithm cleverly exploits this dial. The black box is used to 
produce samples across a range of stepsizes. Most of the samples that we 
ask for will be obtained quickly with relatively large At values. Relatively 
few samples will be generated at the expensive small At levels. In a sense, 
the large At paths cover the low-frequency information so that expensive, 
high-frequency paths are used sparingly. Figure [l] might convince you that 
this idea has some merit. The next section works through the details. 

5 Multilevel Monte Carlo with Euler-Maruyama 

We focus now on the more general case where we wish to approximate the 
expected value of some function of the final time solution, E[h(X(T))]. We 
have in mind the case where X(t) represents an underlying asset, in risk- 
neutral form, and h(-) is the payoff of a corresponding European-style option 
[2D) 23]- For example, h(x) = max(x — E, 0) for a European call option with 
exercise price E and expiry time T. For simplicity, we will consider the scalar 
case, so that m = d = 1 in ([Tj) , but we note that all arguments generalise to 
the case of systems, with the same conclusions. We assume that the payoff 
function h satisfies a global Lipschitz condition; this covers the call and put 
option cases. 

Multilevel Monte Carlo uses a range of different discretization levels. At 
level / we have a stepsize of the form 

Ati = M~ l T, where l — 0,1, 2,..., L. (9) 

Here M > 1 is a fixed quantity whose precise value does not affect the overall 
complexity of the method, in terms of the asymptotic rate as e —)■ 0. For 
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simplicity we may think of M 
choose 


= 2. As the upper limit on the level index we 


loge 1 
log M 


( 10 ) 


In this way, at the coarsest level, l = 0, we have the largest stepsize, At 0 = T, 
covering the whole interval in one step. At the most refined level, l — L, we 
have A t L = 0(e )—from Q, this the stepsize needed by Euler-Maruyama to 
achieve weak error of 0(e). 

With each choice of stepsize, At/, we may apply Euler-Maruyama to the 
SDE (|Tj) and evaluate the payoff function h at the final time. We will let the 
random variable P/ denote this approximation to h(X(T)). Now, from the 
linearity of the expectation operator we have the telescoping sum 


L 

E [P L ] = E[P 0 ] + ^E[P/ - Pi- r]. (11) 

i=i 

In multilevel Monte Carlo, we use the expansion on the right hand side as 
an indirect means to evaluate the left hand side. This may be thought of as 
a recursive application of the control variate technique, which is widely used 
in applied statistics [201 [23U371138]. To estimate E[P 0 ] we form the usual 
sample mean, based on, say, N 0 , paths. This gives 


N 0 




s=l 


( 12 ) 


Generally, for E[P/ — P/_i] 


with / > 0 we will use Ni paths so that 


Y t = 



(13) 


It is vital to point out that P( S and P)_, in (13) come from the same dis¬ 
cretized Brownian path , with different stepsizes At/ and At/_i, respectively. 
Figure [2] illustrates the idea for the case M = 2. In words, at a general 
level l , we compute Ni Brownian paths and, for each path, apply Euler- 
Maruyama twice; once with stepsize At/ and once with stepsize At/_i. (In 
practice, we compute a path at resolution At/ and then combine Brownian 
increments over pairs of steps in order to get a path at resolution At/_j.) 
Having constructed our A/ independent paths for level /, we start afresh at 
level l + 1; none of the earlier information is re-used and new (independent) 
pseudo-random numbers are generated. 




Because of the choice of L in (10) we know from 0 that our estimator 
will have the required 0(e) bias. Now we will see how to choose the values 
of {Ni yf'—Q to achieve the corresponding accuracy in the overall confidence 
interval. 

Considering a general level where l > 0, appealing to the strong con¬ 
vergence behaviour B of Enlcr-Maruyama and our assumption that h is 
globally Lipschitz, we have 

var[fl - fc(A'(T))] = E[(P, - ft(A'(T))) 2 ] - (E[fj - ft(A(T))]) 2 (14) 

< E[(fl - ft(A'(T))) 2 ) (15) 

< constant x E[(Ahv — X(T )) 2 ] (16) 

= 0(At{). (17) 

ft then follows from Minkowski’s Inequality [8] that 

var[P) — Pi-i] = var [P t - h(X(T)) + h(X(T)) - P^] 


< I \/var[P) - h(X(T))} + \/var [h(X(T)) - P z _ x ] 
= O(Ati). 


(18) 


Applying this result in (13) we conclude that Y) has a variance of 0(Ati/Ni ) 
for l > 1. Because all levels are independent, we deduce that 


var 


Y (] 




(=i 


var[Yo] + ^ 0(Ati/Ni). 


i=i 


To balance the variance evenly across levels l = 1, 2,..., L and to control the 
variance at level / = 0, we choose A) = 0(e~ 2 LAti). It then follows that onr 
overall estimator has variance 


0(£ 2 ) + ^0(£ 2 /L)=0(£ 2 ). 


(=1 


In this way, we have achieved the bias and variance required to give a confi¬ 
dence interval of the specified e level of accuracy. 

To quantify how the complexity of this algorithm scales with e, we sum 
the cost of level l from l = 0 to L to give 


J^NiAt; 1 = ^e~ 2 LAt l At; 1 = L 2 e~ 2 . 


1=0 


1=0 


From (10) this expression becomes 0(e 2 (loge) 2 ), as we quoted in section |3 
At this stage, a few remarks are in order: 
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Constructive Upper Bound: In the course of the analysis above, we came 
up with a general-purpose choice for the number of paths at each level, 
{Ni}f : 1 0 - The final complexity count is therefore an upper bound on 
the best possible value. In practice, for a given problem and accuracy 
requirement, we can perform a cheap pre-processing step where appro¬ 
priate variances are estimated and an optimization problem is solved 
in order to give a sequence {Ni}f =0 ] see, for example, [IS] . 


Weak versus Strong: The key inequality (18), which guarantees tight cou¬ 
pling between coarse and fine paths, made use of the strong convergence 
property. For small A t L , both paths are close to the true path, so the 
paths must be close to each other. In this sense, both strong and weak 
error rates are key ingredients in the analysis. We note, however, that 
Giles [13] has also developed estimators that do not rely directly on 
strong convergence. 


Variance and Second Moment In deriving the inequality (17), we dis¬ 
carded the square of the first moment and used the second moment 
as an upper bound for the variance. This may appear to be a very 
crude step, but in our context it does not degrade the final conclusion. 
(In a different, Poisson-driven setting where a multilevel method was 
developed and analysed, the step (14)-(15) is no longer optimal—it is 
beneficial to analyse the variance directly [3].) 


Exploiting Structure: As mentioned above, multilevel Monte Carlo may 
be viewed as a recursive version of the control variate approach. In the 
simplest version of control variates, if we wish to compute E[X], we 
may instead compute E[X — Y] and add E[F], where Y is a suitably 
constructed random variable such that X — Y has small variance and 
E[y] is readily available [201 23 ] • However, the success of this technique 
usually relies on incorporating some extra knowledge of the problem: a 
structure such as symmetry or convexity, or the existence of a “nearby” 
problem that is analytically tractable. In this respect, the multilevel 
Monte Carlo method for SDEs is very different from traditional control 
variates: the analysis is completely general and no special insights are 
needed about the nature of the underlying SDE, other than knowledge 
of the basic weak and strong convergence properties. 


Multilevel versus Multigrid: In [IB], Giles explains that the multigrid 
approach in numerical PDEs was “the inspiration for the author in de¬ 
veloping the MLMC method for SDE path simulation.” There are clear 
similarities between the two: the use of geometrically refined/coarsened 
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grids and the idea that relatively little work needs to be expended on 
the fine grids in order to resolve high frequency components. However, 
it is important to keep in mind that there are also conceptual differ¬ 
ences between the two techniques: multilevel Monte Carlo is distinct, 
and novel. For example, multilevel Monte Carlo does not involve the 
notion of passing information up and down the refinement levels, as is 
done with multigrid V or W cycles. 

Related Earlier Methods: As discussed in pm section 1.3], related ear¬ 
lier work on improving Monte Carlo when samples are generated via 
discretization was performed by Heinrich, see, for example E 2 BE 2 !, and 
Kebaier H devised a two-level approach to path simulation. 

Based on the type of analysis that we summarized above, it is possible to 
state a general theorem about multilevel simulation: 

Theorem 5.1 (Giles; see for example, [ 16 j) Let P denote a random vari¬ 
able, and let Pi denote the corresponding level l numerical approximation. If 
there exist independent estimators Y t based on Ni Monte Carlo samples, and 
positive constants a, P, 7, <7, c 2 , c 3 such that a > \ min(5,7) and 

1. |E[Pj - P] | < Ci2~ al 

2. E[Y 0 ] = E[P 0 ] and E [Yf = E[P, - P_ 4 ] for l > 0 

3. var [Yi\ < c 2 Nf ] 2-^ 1 

4- E[G] < C 3 -/V/ 2 7i , where Ci is the computational complexity ofYi 

then there exists a positive constant c 4 such that for any e < e^ 1 there are 
values L and Ni for which the multilevel estimator 

Y = ± Yl 
1=0 


has a mean-square error with bound 


E [(y - E[P]) 2 ] < e 2 


with a computational complexity C with bound 


E[C] < 


c 4 e 2 , p > 7 

c 4 e~ 2 (log(e)) 2 , ^ = 7 
c 4 e _2_ ^ 7_/3 ^", P < 7. 
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Giles j|13] has also shown how to construct estimators for which /3 > 7 = 
1, by replacing Euler-Maruyama with the more strongly accurate Milstcin 
scheme. For European-style options with Lipschitz payoff functions, this 
makes 0(e -2 ) complexity achievable. From the arguments in section (|3]), it 
is intuitively reasonable that this is the optimal rate. The issue is formalized 
in [35], and optimality is confirmed. 

In section [2] we mentioned that the basic Euler-Maruyama method ([!]) 
may fail to converge in a weak or strong sense on nonlinear SDEs in the 
asymptotic limit At —> 0. A closely related question, of direct relevance 
to this review, is whether the combination of “Euler-Maruyama plus Monte 
Carlo” converges in the e —>■ 0 limit. In [25] 126]. ffutzenthaler and Jentzen 
showed that Euler-Maruyama Monte Carlo can achieve convergence in a P- 
alrnost sure sense in cases where the underlying Euler-Maruyama scheme 
diverges. This can happen when the events causing Euler-Maruyama to di¬ 
verge are so rare that they are extremely unlikely to impact on any of the 
Monte Carlo samples, ffowever, in [28] ffutzenthaler, Jentzen and Platen 
showed that the multilevel Monte Carlo method does not inherit this prop¬ 
erty. They established this result using a counterexample of the form 

dX(t) = -X(t) 5 dt, (19) 


with X(0) having a standard Gaussian distribution, where E[A(f) 2 ] is the 
required moment. Note that the SDE (19) has a zero drift term, so it may 
also be regarded as a random ODE. A modified version of Euler-Maruyama, 
known as a tamed method, was shown in [28] to recover convergence in the 
multilevel setting. 


6 Computational Experiments 

Asymptotic, e —> 0, analysis indicates that multilevel Monte Carlo offers a 
dramatic improvement in computational complexity. Numerous computa¬ 
tional studies have confirmed that this potential can be realised in practice. 

Giles has made MATLAB code available at 

\protect\vrule widthOpt\protect\href{http://people.maths.ox.ac.uk/gilesm/acta/M 

that can be used as the basis for computational experimentation. In Figure [3] 
we show results based on this code. Here, we have an asset model given by 
geometric Brownian motion 

dX{t) = 0.05 X(t)dt + 0.25X(f)dJE(f), X(0) = 100. 
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We consider (a) a European call and (b) a digital call option over [0,T] with 
T — 1 and exercise price 100. So the payoff functions, after discounting for 
interest, are 

h(x) = e _0 05T max(a; — 100, 0) 

for the call option and 


f e °-° 5T 100 when x > 100 
\ 0 when x < 100 


for the digital option. (For those who worry about probability zero events, 
the code defines h(100) = e _0 -° 5T (100 + 0)/2.) The code repeats the Monte 
Carlo simulation for accuracy requests of e = 0.1, 0.05, 0.02, 0.01, 0.005. The 
upper left picture in Figure [3] shows, for the call option, the number of paths 
Ni used at each level l in the multilevel method. We see that for a given e 
more paths are used at the cheaper (small /) levels, and as e is decreased, 
so that more accuracy is required, extra levels are added. The upper right 
picture indicates the corresponding computational cost in terms of run time. 
More precisely, the asterisks (joined by dashed lines) show the cost weighted 
by e 2 as a function of e. We see that this quantity remains approximately 
constant, as predicted by the analysis. The picture also shows the scaled cost 
for an equivalent standard Monte Carlo computation, using a solid linetype. 
We see a much larger cost that appears to grow faster than e -2 . The lower 
pictures in Figure [3] give the same results for the case of the digital option, 
and again the multilevel version is seen to be more efficient than standard 
Monte Carlo. 


7 Follow-on Research 

In this section we summarize some of the key advances that have been made 
since the original multilevel breakthrough D 33 - We focus on work that is 
directly relevant to financial option valuation. The comprehensive overviews 
[IB, EES] can be consulted for further details on these, and other, areas. The 
webpage maintained by Giles at 

\protect\vrule widthOpt\protect\href{http://people.maths.ox.ac,uk/\string~gilesm 

is also an excellent source of up-to-date information. 

7.1 Beyond European Calls and Puts 

A key step in the analysis of section [5] was to show that the coarse and 
refined paths are tightly coupled, in the sense that they produce payoffs 
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whose difference has small variance. The logic behind the analysis may be 
loosely summarized as 

A strong convergence of Euler-Maruyama =>■ 

B coarse and refined paths close to the true path 
C coarse and refined paths close to each other =>- 
D coarse and refined payoffs close to each other. 

The C D step appealed to the global Lipschitz property of h. This is 
valid for European call and put options, where h(x) = max(i — E, 0) and 
h(x) = ma x(E — x,0), respectively. However, the analysis must be refined 
for those European-style options where E[/i(X(T))] is required for functions 
h that violate the global Lipschitz criterion. We may also wish to deal with 
path-dependent options where an expected value operation is applied to a 
functional depending on some or all of the values X(t) for 0 < t < T. 

These more exotic options include problematic classes where, for certain 
SDE paths, the payoff may be very sensitive to small changes. For example, 
with digital options that expire close to the money, a small change in the as¬ 
set path can lead to an 0(1) change in the payoff. Similarly, the payoff from 
a barrier option is very sensitive to those paths that flirt with the barrier. 
In these cases, the logical flow above above must be adapted. Intuitively, 
we should be able to exploit the fact that troublesome paths are the excep¬ 
tion rather than the rule, and hence C D with high probability. In some 
cases this allows us to recover the computational complexity that we saw for 
European calls and puts. In other cases we must accept a slight increase in 
cost. 

The behaviour of multilevel Monte Carlo for Asian, lookback and digital 
options was considered computationally in the original work of Giles [IT]. 
Rigorous analysis to back up these results for barrier, lookback and digital 
options was first given in Id- Further work has been targeted at binary 
options [5], Asian options [2], basket options [T5] , barrier options [10] and 
American options [6]. The use of multilevel Monte Carlo to compute sensi¬ 
tivities with respect to problem parameters, that is, Greeks, was considered 
in [7]. 

7.2 Further Developments 

It is common practice to combine more than one variance reduction tech¬ 
nique. Given that antithetic variables can be effective in option valuation 
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[201 [23], it is natural to consider embedding this approach within the mul¬ 
tilevel framework. Giles and Szpruch HH H9] have shown that this can be 
effective, particularly when Milstein is used for the numerical integration. A 
conditional Monte Carlo approach has also been shown to be fruitful in the 
mutlilevel setting Pi- In a different direction, Rhee and Glynn [36] have pro¬ 
posed an extra level of randomization that produces an unbiased multilevel 
estimator. 

To go beyond the asymptotic 0(e~ 2 ) complexity barrier it is possible 
to move to quasi Monte Carlo, where a low-discrepancy sequence replaces 
a pseudo-random sequence. Giles and Waterhouse have demonstrated 
that a combination of quasi Monte Carlo and multilevel can outperform each 
separate technique. 

Finally, we note that the multilevel methodology has also been extended 
to asset models that are not driven purely by Brownian motion mm- 


8 Discussion 

Our aim in this article was to explain in an accessible manner the key ideas 
behind the multilevel Monte Carlo method. We focussed on the case of 
SDE-based financial option valuation, where Monte Carlo is a widely used 
tool. At the heart of the technique is a very general and widely applicable 
philosophy—a recursive application of control variates that relics on tight 
coupling between simulations at different resolutions. The resulting algo¬ 
rithm is sufficiently simple and effective that it can be implemented straight¬ 
forwardly and used to produce noticable gains in computational efficiency in 
very general circumstances. However, as evidenced by the wealth of current 
research activity, there is also substantial scope for (a) refining the multilevel 
approach in order to exploit problem-specific information and (b) developing 
multilevel methods in many other stochastic simulation scenarios. For these 
reasons we envisage multilevel Monte Carlo evolving into a cornerstone of 
computational finance. 
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Figure 1: Paths based on the Paley- Wiener representation ([8]). As indicated, 
the six plots show the sum truncated after M — 1, 2, 5, 10, 50 and 200 sine 
terms. 
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Figure 2: Illustration of how the estimator Yj in (13) is constructed. Circles 
(joined by straight lines for clarity) show the refined Enler-Marnyama path, 
with stepsize At = 2 ~ l T. Asterisks show the coarser Enler-Marnyama path, 
with stepsize At = 2 ~ l+1 T, computed with the same Brownian increments. 











0 2 4 6 8 10 - 2 10 -' 


l “ v “ t 1 accuracy e 

Figure 3: Output from the multilevel Monte Carlo code made available by 
Giles (see text for web site address). Left hand pictures show the number 
of paths per level at each target accuracy. Right hand pictures show the 
computation time, scaled by e 2 . Upper pictures are for a European call 
option. Upper pictures are for a digital option. 
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